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INTRODUCTION 

There  are  many  users  of  satellite  ephemerides.  Users  of  data  &om  a satellite  such  as 
LANDSAT  need  the  satellite’s  position  to  correlate  the  data  with  a specific  location  on  the 
earth.  The  ephemerides  of  navigation  satellites  are  stored  onboard  and  transmitted  to  users. 
Ground  stations  need  ephemerides  for  antenna  pointing.  In  orbit-prediction  programs,  the 
ephemerides  of  the  planets  and  moon  are  needed  to  compute  their  effects  on  satellite  orbits. 

An  ephemeris  can  be  provided  in  several  ways.  It  can  be  in  tabular  form,  in  which  the 
user  interpolates  between  points  to  obtain  positions.  Orbital  elements  (such  as  NAVSPASUR 
1-  or  5-card  element  sets)  can  be  provided  from  which  the  ephemeris  is  reconstructed  by 
using  a Brouwer  Theory  orbit  generator.  Another  method  is  to  represent  the  ephemeris 
using  polynomial  approximation.  In  this  method,  the  user  is  supplied  with  the  degree  of  the 
approximating  polynomial  and  the  coefficients  needed  to  construct  the  polynomial. 

This  program  began  from  the  desire  to  develop  an  efficient  method  for  storing  and 
processing  satellite  ephemerides  onboard  satellites.  Participation  of  the  senior  author  in  the 
National  Bureau  of  Standards  experiment  of  broadcasting  time  from  geosynchronous  satel- 
lites [1-3]  generated  the  problem  of  loading  astronomical  ephemerides  in  the  memory  of 
microprocessors  and  microcomputers.  The  specific  goals  are 

1.  To  compress  the  ephemerides,  thereby  reducing  the  occupancy  in  core  and  the 
duration  of  loading  into  core, 

2.  To  cover  as  wide  a time  span  as  the  recipient  equipment  will  admit,  in  order  to 
confer  it  with  maximum  autonomy, 

3.  To  guarantee  as  much  accuracy  as  the  real-time  operations  will  require, 

4.  To  make  the  processing  of  the  compressed  ephemerides  as  fast  as  possible. 
Polynomial  approximation  with  Chebyshev  polynomials  can  achieve  all  of  these  goals. 


APPROXIMATION  TECHNIQUES 
Inten>olation  Polynomials 

Polynomial  interpolation  may  help  somewhat  in  compressing  conventionally  tabulated 
almanacs.  Newhall  [4]  reports  how  an  interpolation  in  Chebyshev  polynomials  helped  in 
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making  the  files  of  the  JPL  Development  Ephemeris  Number  96  [ 5]  five  times  shorter  than 
their  homolog  in  the  JPL  Ephemeris  Number  69  [6] . However,  for  an  ephemeris  produced 
by  a numerical  integrator  of  order  n and  step  At,  interpolation  polynomials  are  valid  over 
intervals  equal  to  nAt.  We  are  interested  in  polynomials  extending  over  10  to  40  times  that 
interval.  From  a predictor-corrector  of  order  10  at  the  step  of  500  s applied  to  24-h  satellite, 
the  extraction  of  a polynomial  approximation  valid  for  2,  3,  or  even  5 days  is  proposed. 

Least-Squares  Approximation 

Forsaking  interpolation  polynomials,  we  turn  to  approximation  polynomials.  Let  f(T) 
be  one  of  the  ephemeris  elements.  Define  the  error  y(r)  of  a finite  expansion  by  forming  the 
difference 

y(r)  = f(T)  - ^ CjPy(r) 

0<j<  m 

where  pg,  , . . . , p„  are  prescribed  polynomials  while  the  coefficients  Cp,  c^, . . . , are 
to  be  determined.  Instead  of  forcing  y(r)  to  vanish  at  prescribed  points  (as  an  interp>olation 
requires),  the  average  error  is  obtained  by  integrating  over  the  finite  range  a < r < ^ in 
which  the  function  /(r)  shall  be  represented. 

Ck>nventional  astronomers  usually  define  the  average  error  by  the  definite  integral 

y"  yHr)  dr 

^ •'Of 


and  determine  the  coefficients  Cj  of  the  expansion  by  the  prescription  that  be  as  small  as 
possible.  Least-squares  approximations  always  have  a definite  solution  since  an  essentially 
positive  quantity  always  assumes  a definite  absolute  minimum. 

The  problem  is  greatly  facilitated  if  the  base  polynomials  satisfy  orthogonality  condi- 
tions of  the  kind 


/ 


PjiT)P^iT)  = 0 


Chebyshev  polynomials  are  orthogonal. 


for  /■  ^ k. 


Considerable  effort  has  been  spent  in  the  least-squares  approximation  of  planetary 
ephemerides  in  series  of  Chebyshev  polynomials.  The  primary  example  is  Carpenter’s 
theory  of  the  five  outer  planets  [7] ; following  the  method  outlined  by  Clenshaw  [8,9] , 
Carpenter  develops  a least-squares  approximation  of  the  right-hand  members  of  the  plane- 
tary equations,  then  integrates  formally  the  resulting  Chebyshev  series.  The  approximation 
is  iterative;  starting  with  polynomials  representing  uncoupled  elliptic  Keplerian  motions,  it 
ends  up  with  series  of  several  hundred  terms  representing  the  deviations  of  the  orbits  of 
Jupiter,  Saturn,  Uranus,  Neptune,  and  Pluto  from  Keplerian  motions  in  the  interval  from 
1800  through  2000.  Comparison  with  a numerical  integration  [10]  from  1653  through 
2060  produced  displacements  less  than  0 "01  after  the  third  iteration  and  less  than  0'.'001 
after  the  fourth. 
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Series  of  several  hundred  terms  are  not  practical.  Least-squares  approximations  of 
astronomical  ephemerides  in  Chebyshev  polynomials  are  currently  derived  from  convention- 
aUy  tabulated  ephemerides.  The  orbit  generator  Goddard  Trajectory  Determination  System 
(GTDS)  for  satellites  of  the  earth  or  the  moon  takes  its  almanacs  for  the  sun,  the  moon,  the 
sidereal  time,  and  its  corrections  due  to  precession  and  nutation  from  Chebyshev  series  of 
degrees  which  vary  between  10  and  20.  These  series  are  valid  over  time  intervals  of  20  to  30 
days  and  are  constructed  from  either  the  JPL  files  or  the  Improved  Lunar  Ephemerides  [11]  • 

Techniques  of  this  kind  have  been  in  use  in  producing  ephemerides  which  serve  to 
automatically  drive  antennas  and  signal  transmissions  [12] . They  have  been  transposed  to 
the  problem  of  tracking  specific  features  at  the  moon’s  surface  [13,14] . With  the  advent  of 
programmable  pocket  calculators  available  at  the  drugstore  counter,  the  U.S.  Naval  Observa- 
tory is  now  experimenting  with  a new  type  of  almanac  where  the  major  tables  of  the 
American  Ephemeris  and  Nautical  Almanac  are  compressed  into  standardized  Chebyshev 
series  [15] . 

Chebyshev  Approximations 

The  disadvantage  of  least-squares  approximations  is  usually  the  absence  of  an  estimate 
of  the  error.  The  difference  between  the  function  and  its  polynomial  approximation  in 
least-squares  oscillates  with  nonuniform  magnitude;  normally  minimal  in  the  middle  of  the 
range,  the  amplitude  may  greatly  increase  toward  both  ends  of  the  interval.  A uniform 
approximation  produces  an  error  whose  absolute  maxima  and  minima  are  equal  in  absolute 
value  over  the  entire  interval.  With  a uniform  representation,  there  is  no  need  to  anticipate  a 
“transient”  regime  at  the  beginning  of  the  range  wherein  the  maximum  possible  error  de- 
creases, and  a “degradation”  regime  at  the  end  of  the  range  where  one  should  expect  the 
maximum  error  to  start  increasing  again  in  absolute  value. 

The  topology  subjacent  to  the  Chebyshev  approximation  is  the  metric 

The  purpose  of  the  approximation  algorithm  is  to  find  the  coefficients  Cy  for  which  y* 
would  be  minimum.  The  method  is  based  on  a conjecture  of  Chebyshev  proved  by  Borel 
[16] . If  g(3c)  is  a function  continuous  in  the  closed  interval  (a,  b),  then  there  exists  a unique 
polynomial  of  best  approximation  of  given  degree  n.  The  error  of  this  approximation  reaches 
its  extreme  value  at  n + 2 points  at  least  (there  may  be  more)  with  alternating  signs  at  these 
n + 2 points. 

Stiefel  [17]  has  set  up  an  iterative  method  to  yield  a best  approximation  in  a finite 
number  of  steps.  Assume  that  a sequence  (t,)  of  n + 2 points  has  been  selected  in  the  inter- 
val and  that  a polynomial  series 


<l>iT)  = 2 

0<J<  n 
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has  been  obtained  such  that  the  errors 


y,  = f{r^)  - 0(T,.) 


1 < I < n + 2 


have  the  same  absolute  value,  say  y®.  Consider  the  maximum  y*  over  the  whole  interval.  If 
y * = y® , then  the  polynomial  series  is  the  best  approximation;  otherwise  y*>y^ . Then 
select  one  of  the  points  t at  which  the  maximum  error  is  reached,  exchange  it  with  one  of 
the  points  in  the  reference  set,  and  go  back  to  constructing  a levelled  approximation,  i.e., 
one  for  which  the  errors  at  the  reference  points  are  all  equal  in  absolute  value,  but  alternate 


in  sign. 


Stiefel’s  algorithm  will  confront  the  programmer  with  two  difficulties:  (a)  levelled 
approximations  are  easy  to  construct  [18] , but  only  if  Haar’s  condition  is  satisfied  at  the 
reference  points;  (b)  locating  the  extrema  of  the  error  functions  outside  the  reference  points 
requires  that  one  is  capable  of  determining  accurately  the  first  derivative  f'ir).  If  the  func- 
tion f is  known  in  literal  form,  calculating  f'(T)  usually  presents  no  problem. 

A program  (in  PL/I,  IBM  Level  F)  to  implement  Stiefel’s  iteration  [19,20]  has  been 
developed,  and  it  has  been  applied  to  those  astronomical  ephemerides  which  are  derived 
from  theories,  namely 

1.  the  corrections  due  to  nutation  from  Woolard’s  theory, 

2.  the  coordinates  of  the  sun.  Mercury,  Venus,  and  Mars  from  Newcomb’s  theory  of 
the  inner  planets, 

3.  the  coordinates  of  the  Galilean  satellites  of  Jupiter  from  Sampson’s  tables. 

The  results  are  definitely  encouraging.  For  instance,  polynomials  of  degree  14  are  sufficient 
to  cover  a sidereal  period  (about  230  days)  of  Venus  at  the  accuracy  of  the  American 
Ephemeris.  Similar  results  are  obtained  for  Mars.  The  situation  is  less  favorable  for  Mercury. 
Our  experience  with  the  Galilean  satellites  suggests  that  the  error  which  an  approximation 
of  given  degree  generates  is  determined  for  its  principal  part  by  the  average  eccentricity  of 
its  orbit.  This  correlation  will  be  analyzed  more  closely  later  on. 

A case  where  Haar’s  condition  was  violated  has  not  been  encountered. 

The  precision  of  the  first  and  second  derivatives  is  critical.  A straightforward  derivation 
of  the  semianalytical  expressions  given  by  a theory  is  safe,  but  a numerical  differentiation 
by  finite  differences  is  not  [21] . 


DISCRETE  CHEBYSHEV  APPROXIMATIONS 

A problem  arises  if  the  Chebyshev  approximation  technique  described  in  the  previous 
section  is  applied  to  orbits  of  artificial  satellites  around  the  earth.  These  orbits  are  usually 
given  in  the  form  of  tables  at  equidistant  times  because  they  are  derived  from  a numeric^ 
integration  and  a literal  development.  No  reliance  can  be  made  on  the  Steifel-Remez  algo- 
rithm to  compress  such  ephemerides.  The  construction  proposed  by  Golub  [22]  would  not 
produce  with  enough  accuracy  the  dates  at  which  Remez’s  ripplings  occur. 
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The  following  compromise  was  found  to  be  successful:  in  the  range  -1  < t < +1, 
select  a sequence  ^ and  solve  the  overdetermined  linear  system 

At,)  = {0  < i < m,  n < m) 

0<j<  m 

in  the  sense  of  Chebyshev,  i.e.,  determine  the  coefficients  (c.)  which  minimize  the  error 
estimate  ■' 

**■*'''"  0<><n 

To  accomplish  this,  a very  efficient  algorithri  established  by  Barrodale  [23,24]  can  be  used. 
It  is  basically  the  simplex  algorithm  applied  to  the  dual  of  the  lin^-r  programming  defined 
by  the  minimization  of  the  maximum  error.  At  first  a simplicial  basis  is  built  and  a starter  is 
determined  to  enter  a cycle  of  exchanges  which,  in  a finite  number  of  steps,  leads  eventually 
to  the  optimum  sequence 

Barrodale ’s  program  is  written  in  the  Fortran  dialect  WATFOR.  The  authors  of  this 
report  have  modified  it  to  relax  the  syntactical  restrictions  imposed  by  WATFOR,  to  satisfy 
the  stylistic  rules  adopted  in  structured  programming,  and  to  insert  it  in  production  pro- 
-ams which  compress  by  segments  a conventionally  tabulated  ephemeris  stored  in  a sequen- 
tial file.  The  program  is  available  upon  request  to  the  office  of  the  Space  Systems  Division  at 
the  U.S.  Naval  Research  Laboratory.  Its  input  dispositions  have  been  developed  to  process 
large  ephemerides  files;  this  version  is  available  at  the  Charles  Stark  Draper  Laboratory. 

Illustrations 

As  a test  of  what  discrete  Chebyshev  approximations  can  achieve  in  compressing 
ephemerides,  consider  the  ephemerides  of  the  moon. 

In  principle,  a continuous  Chebyshev  approximation  can  be  used.  The  lunar  ephemeris 
may  be  calculated  from  the  semianalytical  theory  known  as  the  Improved  Lunar  Ephemeris 
(type  j=  2).  But  the  multivariate  Fourier  series  contains  several  hundred  terms,  and  it  is  not 
practical  that  such  a large  series  be  evaluated  repeatedly  not  only  for  the  coordinates  (as 
they  must  be  in  any  case)  but  also  for  their  first  and  second  derivatives  with  respect  to  the 
time. 


Dr.  Van  Flandem  of  the  U.S.  Naval  Observatory  provides  a program  which  optimizes 
the  tabulation  of  a lunar  ephemeris  at  equidistant  times  from  the  amended  series  in  the 
Improved  Lunar  Ephemeris.  The  primary  interest  is  in  testing  how  wide  a range  could  be 
covered  whUe  maintaining  the  precision  retained  in  the  American  Ephemeris  and  Nautical 
Almanac.  The  errors  in  absolute  value  predicted  by  the  algorithm  are  presented  in  Table  1. 
The  range  is  expressed  in  days;  the  deviations  AU,  AV,  Aa,  and  A8  in  longitude,  latitude, 
right  ^cension,  and  declination  are  expressed  in  units  equal  to  10'®  radian,  while  the  error 
Ar  is  in  units  of  10'®  times  the  earth  radius.  The  table  shows  that  a full  period  of  the 
moon  may  be  covered  by  Chebyshev  polynomials  of  only  degree  24  at  a precision  of  0'.'035 
in  the  angles  and  of  16  m in  the  distance  (comparable  to  the  accuracy  maintained  by  the 
American  Ephemeris  and  Nautical  Almanac). 
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Range  Degree 


Table  I — Compression  of  Lunar  Ephemeris 


Ar 

Aa 

A6 

28897 

4840 

5674 

529 

1706 

635 

437 

401 

157 

323 

40 

7 

316 

53 

11 

264 

46 

11 

226 

43 

10 

198 

40 

10 

131 

37 

9 

102 

35 

9 

90 

32 

8 

576 

630 

603 

164 

823 

806 

3133 

6746 

2752 

2213 

1901 

357 

1988 

377 

132 

2453 

5425 

1337 

2321 

5196 

699 

2791 

525 

136 

2471 

171 

35 

1274 

25 

6 

1557 

34 

8 

1457 

32 

7 

3304 

1317 

1019 

291736 

184845 

141263 

4493 

2998 

2773 

2170 

809 

109 

1549  1 

94  1 

36 

A similar  analysis  of  possible  ranges  and  degrees  was  made  for  the  satellite  SMS-B  of 
NOAA.  The  orbit  was  produced  by  GTDS,  the  geocentric  coordinates  being  filed  at  the  rate 
of  one  state  vector  every  15  min.  A shorter  filing  step  would  have  improved  the  accuracy  of 
the  interpolation  in  the  data  set.  The  best  residuals  predicted  by  Barrodale’s  algorithm 
turned  out  to  be  comparable  to  the  interpolation  errors.  Table  2 suggests  the  type  of  per- 
formance that  one  may  expect  from  the  compressed  ephemerides.  At  the  National  Bureau 
of  Standards,  in  the  present  stage  of  development,  accuracy  requirements  are  moderate 
without  being  acute,  but  it  is  critical  to  secure  as  long  an  autonomy  of  the  equipment 
installed  at  Wallops  Island  as  is  compatible  with  a time  broadcast  to  better  than  10  fxs  to 
general  users.  In  that  context,  it  is  quite  encouraging  to  read  in  Table  2 that  several  days  of 
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Table  2 — Compression  of  Ephemerides 
for  the  24-h  Satellite  SMS-B 


Range 

Degree 

Ar 

— 

AU 

AV 

0‘*.25 

4 

575 

43 

310 

8 

3 

3 

4 

12 

2 

2 

1 

16 

0 

0 

0 

0‘*.50 

8 

228 

9 

5 

12 

5 

3 

4 

16 

4 

3 

3 

20 

2 

2 

2 

0‘^.75 

8 

860 

64 

43 

12 

9 

6 

6 

16 

5 

4 

4 

20 

4 

3 

4 

12 

83 

11 

5 

16 

5 

5 

7 

20 

5 

4 

5 

24 

5 

4 

4 

2^. 

16 

3165 

401 

19 

20 

347 

56 

9 

24 

55 

12 

6 

28 

15 

5 

6 

20 

11676 

2019 

110 

24 

2424 

502 

37 

28 

617 

124 

13 

32 

226 

31 

7 

4^. 

24 

36057 

6433 

389 

28 

5984 

1068 

67 

32 

2974 

529 

39 

36 

968 

201 

18 

28 

67211 

13197 

798 

32 

16257 

2718 

186 

36 

4653 

780 

55 

40 

3345 

585 

42 

DEPRIT,  PICKARD  AND  POPLARCHEK 


Iephemerides  may  be  compressed  in  Chebyshev  series  of  reasonable  degree.  There,  the 

residuals  between  data  from  integration  and  values  from  minmax  approximations  are  in  cm 
for  the  distance  r and  in  degrees  for  longitude  U and  latitude  V. 

In  Fig.  i,  the  error  in  distance  is  plotted  for  an  approximation  of  degree  20  over  two 
I days.  The  curve  is  typical  of  a discrete  Chebyshev  approximation.  The  rippling  character  is 

well  in  evidence;  maxima  and  minima  are  not  quite  at  the  same  height,  which  shows  that  the 
Chebyshev  series  narrowly  misses  being  the  best  appro  limation  of  degree  20  for  the  distance 
over  two  days. 
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Selection  of  the  Reference  Points 

High-frequency  oscillations  of  the  truncated  series  around  the  true  function  are  always 
present.  With  a least-squares  approximation  they  are  very  noticeable  and  interfere  most 
adversely  with  an  efficient  polynomial  synthesis.  The  error  plot  in  Corio  [25]  is  a striking 
illustration  of  Gibbs  phenomenon;  for  a 24-h  satellite,  an  error  of  less  than  1 m around 
noon  causes  wlix.'  oscillations  of  10  km  on  both  sides  of  midnight.  Corio ’s  figure  is  in  sharp 
contrast  with  the  error  curve  in  Fig.  1. 


Fejer’s  rule  of  taking  the  arithmetic  mean 


where 


1 

n + 1 


0<j<n 


Sy  = S (0<y<n) 

o^jifcy 


pliminatpg  completely  the  Gibbs  phenomenon  [26] . Yet  it  must  be  said  that  the  method 
approaches  the  limit  value  very  slowly.  The  method  of  smoothing  the  coefficients  of  the 
approximation  [27,28]  is  a somewhat  simpler  process;  it  does  not  eliminate  the  oscilla- 
tions but  cuts  dovm  their  amplitudes. 


On  the  whole,  however,  much  can  be  accomplished  toward  damping  the  Gibbs  oscilla- 
tions by  securing  a nonuniform  distribution  of  the  approximation  points  which  crowds  the 
points  near  the  two  end  points  t = ±1  of  the  range.  The  effectiveness  of  this  method  is 
illustrated  in  Fig.  2,  (a)  and  (b).  Plotted  here  are  error  curves  corresponding  to  polynomial 
approximations  of  a 12-h  satellite  of  eccentricity  0.01  and  inclination  63.4°  over  two  peri- 
ods. Twenty-five  points  were  used  to  generate  the  approximations  of  Fig.  2(a)  and  sixty 
points  were  used  for  those  of  Fig.  2(b).  The  polynomial  approximations  corresponding  to 
the  three  plots  of  each  figure  were  generated  as  follows: 

Plot  1 : Reference  points  were  located  at  the  zeros  of  the  Chebyshev  polynomial  of 
order  25  (Fig.  2(a))  and  60  (Fig.  2(b)).  The  Barrodale  algorithm  was  used. 


Plot  2:  Reference  points  were  uniformly  distributed.  The  Barrodale  algorithm  was 

used. 


Plot  3:  Reference  points  were  located  as  in  Plot  1,  but  least-squares  fitting  was  used. 


The  contrast  in  results  between  a uniform  and  a nonunifoim  distribution  of  reference 
points  is  most  striking  when  only  25  points  are  used,  but  the  effect  is  still  s^ificant  for  60 
points.  In  Fig.  2(a),  the  accuracy  of  the  approximation  corresponding  to  a uniform  distribu- 
tion deteriorates  steadily  and  drastically  fiom  the  midpoint  to  the  end  points  of  the  interval. 
The  case  for  the  uniform  distribution  is  not  nearly  so  bad  if  60  points  are  used,  but  errors  at 
the  very  end  points  still  extend  by  a factor  of  10  above  those  corresponding  to  nonuni- 
formly  distributed  points.  The  third  plot  in  each  figure  is  included  for  comparison  between 
the  least-squares  fitting  and  the  Barrodale  algorithm.  The  results  are  essentiaUy  the  same  for 
these  two  methods.  The  error  curve  of  a least-squares  fit  for  uniformly  distributed  points  is 
very  similar  to  the  second  plot  and  has  not  been  included  in  Fig.  2. 
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TIME  IN  HOURS 

(a)  Twenty-five  reference  points 


TIME  IN  HOURS 
(b)  Sixty  reference  points 
Fig.  2 — Error  curves  for  the  24th-degree  fit 


Only  when  a large  number  of  reference  points  is  used  can  the  results  of  a uniform  dis- 
tribution match  the  results  of  a nonuniform  distribution.  Figure  3 shows  how  well  the  Bar- 
rodale  maxunum  fitting  error  estimates  the  true  maximum  error  over  the  entire  interval  as  a 
function  of  the  number  of  reference  points  used.  The  solid  curves  indicate  true  errors  and 
the  dashed  curves  indicate  errors  given  by  the  Barrodale  routine.  Curves  1 correspond  to  the 
case  where  the  reference  points  are  taken  at  the  zeros  of  the  Chebyshev  polynomial  of  order 
given  by  the  x-axis.  Curves  2 correspond  to  the  case  where  the  reference  points  are  uniform- 
, The  12-h  satellite  orbit  used  previously  was  also  used  in  generating  these 

plots.  For  these  cases,  it  takes  60  to  100  uniformly  distributed  points  to  match  the  accuracy 
pven  by  25  nonuniformly  generated  points.  Furthermore,  although  the  predicted  error  is 
less  for  uniformly  distributed  points  than  for  nonuniformly  distributed  points,  the  true 

uniform  distribution,  even  if  a large  number  of  points  is'  used  Thus 
while  the  true  error  is  underestimated  in  both  cases,  reliable  prediction  is  at  least  possible  ’ 

even  with  only  25  reference  points  when  they  are  concentrated  near  the  end  poinU  of  the 
intervaJ. 

1 * "owding  of  the  reference  points  toward  the  end  points,  the  error  oscil- 

lates with  the  same  order  of  magnitude  throughout  the  range.  The  error  profile  for  a least- 
^uares  approximation  is  now  similar  to  that  of  a discrete  Chebyshev  approximation.  The 
latter,  however,  would  still  be  preferable  since  it  produces  an  estimate  of  the  least  maximum 
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NUMBER  OF  REFERENCE  POINTS 

Fig.  3 — Comparison  of  polynomial  approximations  using  uniform  and  nonuniform 
reference-point  distributions  for  several  numbers  of  reference  points 


error  that  may  be  reached  for  a given  degree.  This  specific  information  is  critical  in  the 
design  of  a microprocessor  or  in  the  programming  of  a minicomputer  to  process  compressed 
ephemerides. 

Ck>nvergence  and  Eccentricity 

To  get  an  idea  of  how  well  a Chebyshev  series  can  approximate  orbit  functions,  the 
Barrodale  algorithm  has  been  used  to  generate  approximating  polynomials  for  elliptic  orbits 
of  different  eccentricities.  The  following  discussion  and  accompanying  tables  give  a rough 
estimate  of  what  order  is  needed  to  approximate  a given  Keplerian  element  to  within  a 
specified  error. 
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Tables  3-6  show  results  from  a 12-h  orbit  of  various  eccentricities.  The  inclination  was 
63.4°  and  the  argument  of  perigee  was  on  the  equator.  The  approximating  polynomials  were 
generated  by  using  60  reference  points  over  the  interval  of  consideration;  the  points  were 
chosen  at  the  zeros  of  the  Chebyshev  polynomial  of  degree  60.  The  maximum  error  in  fitting 
these  60  points  was  given  by  the  Barrodale  algorithm,  and  this  error  was  used  initially  to 
determine  the  lowest  degree  needed  to  fit  the  entire  interval  to  within  the  specified  error. 
The  resulting  approximating  polynomial  was  then  compared  to  the  actual  function  at  500 
points  over  the  interval  to  determine  the  true  maximum  error.  The  maximum  error  given  by 
the  Barrodale  algorithm  is  always  optimistic  to  some  degree,  and  in  some  cases  this  causes 
the  minimum  degree  to  be  underestimated.  Experience  so  far  indicates  that  for  small  to 
moderate  eccentricities  the  maximum  error  at  the  60  points  is  optimistic  by  no  more  than 
about  10  to  20  percent.  Thus,  only  when  this  error  falls  very  close  to  the  maximum  error 
specified  is  it  necessary  to  increase  the  order  to  ensure  the  desired  accuracy  over  the  entire 
interval.  It  is  stressed,  however,  that  in  general,  this  holds  only  when  the  reference  points  are 
chosen  with  a concentration  near  the  end  points  of  the  interval.  If  a uniform  distribution  of 
points  is  used,  it  takes  considerably  more  points  to  achieve  the  same  accuracy.  In  the  tables, 
an  asterisk  indicates  that  the  degree  needed  is  greater  than  59. 

Tables  3-6  show  results  for  the  radius  component  and  for  the  x-component  in  an  Earth- 
fixed  cartesian  coordinate  system.  Results  for  the  y-  and  z-components  are  analogous  to 
those  for  x.  Problems  were  encountered  in  fitting  latitude  and  longitude,  and  these  are  dis- 
cussed separately. 

Experience  suggests  the  following  conclusions  in  regard  to  fitting  the  cartesian  compo- 
nents of  elliptic  orbits.  Large  eccentricity  causes  the  greatest  difficulty  in  fitting  orbit  func- 
tions. The  degree  needed  to  achieve  a specified  accuracy  becomes  sensitive  to  eccentricity 
when  it  exceeds  0.01  and  extremely  sensitive  when  it  gets  above  0.1.  Not  only  does  the 
necessary  degree  become  large,  but  the  accuracy  improves  very  slowly  with  increasing  degree. 
As  Tables  3 and  4 show,  the  increase  in  degree  needed  to  improve  the  accuracy  by  10'^  is 
only  4 for  e = 0.001,  but  jumps  to  20  for  e = 0.5  and  to  40  for  e = 0.75.  Extending  the 
time  interval  of  consideration  also  has  a great  effect  on  the  degree  needed  to  fit,  but  this  is 
not  a problem  inherent  to  the  particular  coordinate  of  interest.  The  pertinent  question  in 
this  case  is  whether  or  not  it  is  advisable  to  fit  one  large-degree  polynomial  rather  than  two 
smaller  degree  polynomials  to  two  or  more  periods  of  an  orbit.  Tables  5 and  6 show  that  for 
small  eccentricities,  two  periods  can  be  fitted  with  less  than  double  the  degree  needed  to  fit 
one  period.  For  e = 0 ;5,  however,  more  than  triple  the  degree  is  needed  for  two  periods. 
The  general  conclusion  appears  to  be  that  the  cartesian  coordinates  of  orbits  can  be  fitted 
very  well  if  the  eccentricity  is  small,  but  with  difficulty  when  the  eccentricity  is  large. 

The  discussion  above  still  holds  when  one  tries  to  fit  either  latitude  or  longitude,  except 
that  the  inclination  has  a very  definite  effect  on  the  results.  In  fact,  the  longitude  for  a 
63.4°  inclined  orbit  is  already  very  difficult  to  fit,  so  that  results  in  Table  7 are  for  a 10° 
inclined  orbit.  This  sensitivity  to  inclination  is  attributed  to  the  very  rapid  change  of  longi- 
tude at  high  latitudes.  Aside  from  the  sensitivity  to  inclination,  fitting  the  longitude  is  gen- 
erally more  difficult  than  fitting  the  cartesian  components  and,  therefore,  results  are  shown 
for  only  one  period.  Fitting  the  longitude  for  orbits  of  inclination  greater  than  45°,  even  for 
one  period,  is  deemed  impractical.  When  it  is  desired  to  fit  the  longitude  of  an  inclined  orbit, 
the  procedure  should  be  to  first  rotate  the  coordinate  system  so  that  the  x-y  plane  coincides 
with  the  orbital  plane. 


12 


NRL  REPORT  8280 

Table  3 — Minimum  Degree  Necessary  to  Fit  Radius  Components  of 
Elliptic  Orbits  Over  One  Period 


0.001+ 


Max.  Error 


10  km  0 4 

1 km  0 4 

100  m 0 6 

10  m 0 8 

1 m 0 8 

Indicates  that  the  degree  needed  is  greater  than  59. 
Eccentricity 


Degree 


0.01+ 


0.1+ 

0.5+ 

0.75+ 

6 

12 

28 

8 

18 

30 

12 

24 

42 

12 

26 

48 

16 

34 

Sc 

Table  4 — Minimum  Degree  Necessary  to  Fit  X-Component  of 
Elliptic  Orbits  Over  One  Period 


Degree 


0.01+ 


Max.  Error 


10  km 
1 km 
100  m 
10  m 
Im 


Indicates  that  the  degree  needed  is  greater  than  59. 
Eccentricity 


Table  5 — Minimum  Degree  Necessary  to  Fit  Radius  Component  of 
Elliptic  Orbits  Over  Two  Periods 


Max.  Error 

0.0+ 

0.001+ 

0.01+ 

0.1+ 

0.5+ 

10  km 

0 

6 

8 

16 

59 

1 km 

0 

8 

12 

22 

100  m 

0 

10 

14 

28 

♦ 

10  m 

0 

12 

18 

36 

« 

1 m 

0 

14 

22 

42 

* 

Indicates  that  the  degree  needed  is  greater  than  59. 
Eccentricity 
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Table  6 — Minimum  Degree  Necessary  to  Fit  X-Component  of 
Elliptic  Orbits  Over  Two  Periods 


Max.  Error 

Degree 

0.0+ 

0.001+ 

0.01+ 

0.1+ 

0.5+ 

0.75+ 

10  km 

16 

16 

18 

24 

* 

* 

1 km 

18 

18 

20 

32 

* 

100  m 

20 

20 

24 

38 

* 

* 

10  m 

22 

24 

28 

45 

* 

* 

1 m 

22 

26 

30 

52 

* 

* 

* Indicates  that  the  degree  is  greater  than  59. 
^Eccentricity 


Table  7 — Minimum  Degree  Necessary  to  Fit  the  Longitude 
of  Elliptic  Orbits  Over  One  Period 


Max.  Error 
(Radians 

X 10'®) 

Degree 

Max.  Error 
(Sec.  of 
Arc) 

0.0+ 

0.001+ 

0.01+ 

0.1+ 

1000.0 

7 

5 

7 

7 

n 

206.0 

100.0 

9 

9 

9 

9 

20.6 

10.0 

11 

13 

13 

13 

23 

45 

2.06 

1.0 

15 

15 

15 

15 

29 

55 

0.206 

0.1 

17 

17 

17 

17 

35 

♦ 

0.0206 

* Indicates  that  the  degree  needed  is  greater  than  59. 
^Eccentricity 


DERIVATIVE  OF  A CHEBYSHEV  SERIES 

Clenshaw  has  proposed  a recursive  algorithm  to  evaluate  expressions  of  the  type 

m = 5]  c.TfT) 

0</<  n 

in  a straightforward  manner,  i.e.  without  expanding  the  Chebyshev  polynomials.  It  is  defi- 
nitely of  interest  to  supplement  the  procedure  with  one  to  calculate  /"(t)  immediately  from 
f(T)  without  having  to  execute  and  expand  the  derivatives.  This  means  that  in  cases  where 
the  time  rates  of  the  ephemeris  coordinates  would  be  of  use,  one  would  not  have  to  generate 
and  transmit  polynomial  approximations  for  the  derivatives.  The  procedure  is  based  on  the 
following  lemma  [29] : 
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Let  (Py),  {Uj),  and  (Pj)  (for  0 < j)  be  three  sequences  of  functions  such  that 
P + %Po  = 0 

Pn*2  + “n  + lPn  + l + ^nPn  =0  ^ 0)- 


Given  the  sum 


f(r)  = ^ CjPj(t), 

0<J<  n 

construct  the  sequence  (&y)o<;<  n + 2 recurrence  starting  with 


^n*2  ^n*l  ~ ® 


and  looping  through 


bj  + 


for  ;■  decreasing  from  n to  0.  Then  the  value  of  f at  t is  the  sum 
f(T)  = (6q  + o(q(t)6j)Pq(t)  + 6iPi(t). 

When  pj  is  the  Chebyshev  polynomial  of  degree  j,  Clenshaw’s  fundamental  recursive 
identities  are  satisfied  by 


«n  = -T. 


otj  = -2t 


for;  > 1, 


Pj  = l 


for;  > 0. 


Therefore,  the  series  (bj)  is  given  by  the  recursive  equalities 
*’n+2  ’ *n  + l “ ® 


^ - ft. ^2 


for  n > ;■  > 1, 


and  it  yields  the  final  result 


f(T)  = b^T  - ftj  + Cq. 

In  order  to  apply  Clenshaw’s  lemma  to  the  evaluation  of  f'ir),  introduce  the  variable  d 
such  that 

7 = cos  6 
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and  recall  that,  in  terms  of  0 , 


T.  = COS/0. 

Now  putting 

U.(j)  = sin  (/ ■■ 

it  is  found  that 

T'fr)  = /•(/.., 

and  that 

A/)=  X) 

0</<  n 

From  the  obvious  identity 

U-iT)  = 2tUj_ 

U > 0). 


U > 1) 


it  is  found  that  the  functions  U.  verify  the  conditions  of  Clenshaw’s  lemma  for 
Oij  = -2t  and  = 1 for  j > 0. 

Hence  the  sequence 

= 0. 

bj  = 2Tbj_^  — +0  + l)Cy 

leads  to  the  evaluation 

f\T)  = 6o. 

Note  that  the  above  algorithm  is  shorter  and  simpler  than  the  one  proposed  by  Broucke  [30] . 

In  both  procedures,  the  sequence  (6  ) is  not  properly  a one-dimensional  array  (one 
does  not  need  to  save  the  intermediate  values  (6.))  but  a three-ceU  push-dovm  stack.  Thus 
Clenshaw  s procedures  are  particularly  suited  to  microprocessors  or  microcomputers  de- 
signed around  the  chip  of  a pocket  calculator. 

CONCLUSIONS 

The  problem  of  transmitting  and  processing  ephemerides  in  real  time  has  troubled 
astronomere  now  for  some  time.  It  is  generally  thought  that  approximations  in  Chebyshev 
polynomials  are  the  solution  to  the  problem. 
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The  shorter  the  range  of  approximation,  the  lower  the  error.  But  if  the  degree  is  al- 
lowed to  take  somewhat  large  (in  excess  of  10)  values,  the  range  may  easily  extend  over 
several  orbital  periods.  This  contradicts  the  conclusion  drawn  by  VanDierendonck  [31]  that 
approximations  by  polynomials  cannot  cover  more  than  a full  period  of  the  orbit. 

Proper  crowding  of  the  data  points  at  both  ends  of  the  range  damps  the  Gibbs  oscilla- 
tions at  the  extremities  to  the  extent  that  the  error  curve  of  a least-squares  fit  tends  to 
acquire  the  uniform  rippling  character  of  a Chebyshev  (also  called  minmax  or  “best”)  approxi- 
mation. Nevertheless,  discrete  Chebyshev  approximations  are  favored  because  they  rely  on 
global  error  to  iterate  the  calculation  of  the  coefficients  and  thus  produce  directly  an  esti- 
mate of  the  maximum  error  over  the  range.  On  the  whole,  they  produce  shorter  series,  i.e., 
without  tails;  the  range  of  the  coefficients  is  narrower,  and  this  facilitates  the  task  of  eco- 
nomically sizing  the  bit  length  to  stack  the  coefficients  in  the  transmission  message.  On  this 
point  though,  the  reader  should  be  cautioned  against  believing  that  the  authors  advocate  a 
particular  set  of  elements  as  being  the  most  economical  one.  It  so  happens  that,  either  for  the 
24-h  satellites  or  for  the  moon,  the  context  in  which  the  series  would  be  used  requires 
spherical  coordinates.  It  is  likely  that  orbital  elements  rather  than  coordinates  would  be 
more  economical:  the  slow-varying  elements  would  ^ve  rise  to  polynomial  approximations 
of  very  low  degree,  and  only  the  fast-varying  element  would  reach  the  kind  of  high  degrees 
mentioned  previously  in  Tables  1 and  2.  For  orbits  at  low  inclinations  with  small  eccen- 
tricities, the  ideal  coordinates  [32]  or  a variant  thereof  called  the  equinoctial  elements  is 
suggested  [ 33  ] . 
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